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A  continuum  description  of  nonlinear  elasticity, 
slip  and  twinning,  with  application  to  sapphire 

By  J.  D.  Clayton* * 
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A  model  is  developed  for  elasticity,  plasticity  and  twinning  in  anisotropic  single  crystals 
subjected  to  large  deformations.  Dislocation  glide  and  deformation  twinning  are 
dissipative,  while  energy  storage  mechanisms  associated  with  dislocation  lines  and 
twin  boundaries  are  described  via  scalar  internal  state  variables.  Concepts  from 
continuum  crystal  plasticity  are  invoked,  with  shearing  rates  on  discrete  glide  and 
twinning  systems  modelled  explicitly.  The  model  describes  aspects  of  thermomechanical 
behaviour  of  single  crystals  of  alumina  over  a  range  of  loading  conditions.  Resolved  shear 
stresses  necessary  for  glide  or  twin  nucleation  at  low  to  moderate  temperatures  are 
estimated  from  nonlinear  elastic  calculations,  theoretical  considerations  of  Peierls 
barriers  and  stacking  fault  energies  and  observations  from  shock  physics  experiments. 
These  estimates  are  combined  with  the  existing  data  from  high-temperature  experiments 
to  provide  initial  yield  conditions  spanning  a  wide  range  of  temperatures.  The  model 
reflects  hardening  of  glide  and  twin  systems  from  dislocations  accumulated  during  basal 
slip.  Residual  elastic  volume  changes,  predicted  from  nonlinear  elastic  considerations  and 
approximated  dislocation  line  energies,  are  positive  and  proportional  to  the  dislocation 
line  density.  While  the  model  suggests  that  generation  of  very  large  dislocation  densities 
could  influence  the  pressure-volume  response,  volume  increases  from  defects  are 
predicted  to  be  small  in  crystals  deformed  via  basal  glide  on  a  single  system. 
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1.  Introduction 


When  subjected  to  large  deformations  or  mechanical  stresses,  crystalline  solids 
may  deform  by  a  number  of  mechanisms.  In  ductile  crystals,  where  the 
dislocation  mobility  is  sufficiently  high  in  an  adequate  number  of  directions, 
dislocation  glide  is  often  the  predominant  accommodation  mode  for  deviatoric 
deformations  at  stresses  above  the  elastic  limit.  However,  under  different 
circumstances,  inelastic  strains  may  be  accommodated  by  other  mechanisms 
such  as  deformation  twinning  or  fracture. 
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Deformation  twinning  is  thought  to  be  preferred  over  slip  in  cases  where 
resistances  to  dislocation  glide  are  very  large  in  certain  directions  (Hirth  &  Lothe 
1982),  often  in  crystals  of  non-cubic  symmetry.  Furthermore,  twinning  may  be 
preferable  to  slip  at  low  temperatures  in  cubic  crystals,  especially  those  with  low 
stacking  fault  energy  (Bernstein  &  Tadrnor  2004).  As  reviewed  by  Christian  & 
Mahajan  (1995),  a  number  of  theories  exist  that  describe  twin  nucleation  and 
propagation.  Here,  deformation  twinning  is  associated  with  irreversible  shape 
deformation  often  attributed  to  the  motion  of  partial  dislocations  and 
corresponding  formation  of  stacking  faults  (Bilby  &  Crocker  1965;  Christian  & 
Mahajan  1995),  although  twinning  in  a  more  general  sense  may  encompass  a 
variety  of  energy  invariant  transformations  of  the  lattice  (Bhattacharya  1991). 
Although  the  general  theory  developed  in  this  paper  applies  to  any  crystalline 
solid  deformable  via  some  combination  of  elasticity,  slip  and  twinning,  the 
present  application  of  the  theory  focuses  on  ceramic  crystals.  As  defined  here,  a 
ceramic  crystal  is  an  inorganic,  non-metallic  and  ordered  solid.  Electronic 
structures  of  ceramics  feature  ionic  and/or  covalent  bonds  as  opposed  to  metallic 
bonds  found  in  ductile  metals.  The  nature  of  interatomic  forces  in  ionic  and 
covalently  bonded  solids  tends  to  correlate  with  a  high  Peierls  barrier  that 
inhibits  dislocation  motion  at  low  temperatures  (Friedel  1964). 

Single  crystals  of  the  ceramic  alumina  (A1203)  are  known  as  corundum.  The 
stable  phase  at  room  temperature  and  pressure  is  a-corundum,  and  although 
other  phases  exist  (Holm  et  al.  1999),  the  term  corundum  hereafter  refers  to  the 
a-phase.  The  crystallography  of  corundum  can  be  described  via  either 
rhombohedral  or  hexagonal  notation  (Kronberg  1957).  Corundum  is  centro- 
symmetric  and  hence  does  not  exhibit  the  piezoelectric  effect.  Corundum  that  is 
red  in  colour,  from  elemental  traces  of  chromium,  is  known  as  ruby,  with  the 
chromium  content  providing  an  increase  in  the  tensile  yield  stress  at  high 
temperatures  (Klassen-Neklyudova  et  al.  1970).  Corundum  of  all  other  colours  is 
called  sapphire.  Sapphire  is  very  hard,  with  a  value  of  9  on  Mohr’s  scale,  and 
exhibits  an  extremely  high  Hugoniot  elastic  limit  (HEL),  with  values  in  excess  of 
20  GPa  for  single  crystals  of  certain  orientations  (Graham  Sz  Brooks  1971). 
Alumina’s  high  hardness,  high  compressive  strength  and  low  mass  density 
relative  to  steels  render  the  ceramic  a  useful  material  for  many  engineering 
devices  and  systems.  Single  crystals  of  adequate  purity  are  transparent,  with 
widespread  uses  in  optics  and  electronics.  Experimental  data  quantifying  the 
mechanical  response  of  sapphire  single  crystals  of  various  orientations  and 
spanning  a  wide  range  of  temperatures  and  strain  rates  are  available  (Graham  & 
Brooks  1971;  Tyrniak  &;  Gerberich  2007).  While  various  kinetic  relationships  for 
slip  (Lagerlof  et  al.  1994;  Rodriguez  et  al.  2008)  and  twinning  (Scott  &  Orr  1983) 
have  been  postulated  to  describe  individual  experiments,  a  more  general  model  is 
needed  to  collectively  explain  the  material’s  behaviour  over  a  broad  range  of 
loading  conditions  encountered  in  practical  applications. 

The  remainder  of  this  paper  is  organized  as  follows.  In  §2,  physical 
descriptions  of  elasticity,  plasticity  and  twinning  in  crystalline  solids  are  given, 
in  order  to  clarify  distinctions  among  the  three  mechanisms  and  provide 
sufficient  basis  for  the  theoretical  framework  of  §3.  This  framework  is  general 
enough  to  apply  to  any  crystalline  solid  that  undergoes  large  deformations  via 
elasticity,  plasticity  and/or  mechanical  twinning.  In  §4,  the  model  is  specialized 
to  describe  the  behaviour  of  sapphire.  Crystal  structure,  physical  properties, 
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yield  mechanisms,  strain  hardening,  defect  accumulation  and  irreversible  volume 
changes  associated  with  lattice  defects  are  considered  in  turn. 

The  following  notational  scheme  is  used.  Scalars  and  individual  components  of 
vectors  and  tensors  are  written  in  italic  font,  while  vectors  and  tensors  are 
written  in  bold-italic  font.  Einstein’s  summation  convention  applies  for  repeated 
indices.  The  •  symbol  denotes  the  scalar  product  (a-b  =  aaba  =  a1b1  +  a2b2  + 
a3 b3),  while  <g>  indicates  the  outer  product  ((a<E>b)ab  =  aabh).  Juxtaposition 
of  second-rank  tensors  implies  summation  over  one  set  of  adjacent  indices 
(AB— AabBi,c).  Summation  over  two  sets  of  adjacent  indices  is  denoted  by  : 
(A  :  B=AabBab).  Indices  m  parentheses  are  symmetric  (2./1^q^  t—  A^^), 

indices  in  braces  are  antisymmetric  {2A^  =  Aab  —  Aba).  Superposed  ■,  T  and  —  1 
denote  a  material  time  derivative,  transposition  and  matrix  inversion. 
Subscripted  commas  denote  partial  coordinate  differentiation. 

2.  Background:  deformation  mechanisms 

(а)  Elasticity 

A  crystal  is  said  to  deform  elastically  in  the  absence  of  generation  or  motion  of 
defects.  At  the  atomic  scale,  elastic  deformation  alters  relative  distances  and/or 
orientations  among  neighbouring  atoms  within  each  crystallographic  unit  cell 
(Born  &  Huang  1954).  Any  resulting  changes  in  interatomic  forces  produce 
mechanical  stresses  when  the  crystal  is  viewed  as  a  continuous  solid.  Removal  of 
mechanical  stresses  restores  the  original  lattice  without  dissipation  of  energy; 
hence,  elastic  deformation  is  said  to  be  thermodynamically  reversible.  Here,  elastic 
deformation  also  includes  changes  in  interatomic  bond  vectors  induced  by  changes 
in  temperature.  Increases  in  thermal  energy,  i.e.  local  atomic  vibrations,  may 
correlate  with  the  expansion  of  the  lattice  in  the  absence  of  mechanical  stresses. 

(б)  Plasticity 

Plastic  deformation  as  defined  here  occurs  via  glide  of  dislocations  of  edge, 
screw  and/or  mixed  character,  including  loops,  and  encompassing  cross-slip  but 
not  climb.  As  dislocations  travel  through  a  region  of  the  lattice,  the  shape  of  the 
material  will  change,  but  the  interatomic  bond  vectors  remain  the  same,  so  long 
as  no  defects  are  left  behind  within  the  region.  Mechanical  stresses  are 
conventionally  required  to  enact  the  net  glide  of  dislocations  (apart  from 
random  thermal  fluctuations),  for  example  resolved  shear  stresses  must  exceed 
the  Peierls  barrier  in  the  context  of  lattice  statics  or  Schmid’s  law  in  the  context 
of  continuum  slip  (Hirth  Sz  Lothe  1982).  Plasticity  is  thermomechanically 
irreversible,  since  the  reference  shape  of  the  material  is  not  recovered  upon 
removal  of  mechanical  stresses,  and  heat  is  dissipated  by  moving  dislocations  as  a 
result  of  lattice  friction,  phonon  drag  and  other  mechanisms  (Kocks  et  al.  1975; 
Gilman  1979). 


(c)  Twinning 

Deformation  twinning  results  in  two  connected  regions  in  the  lattice  separated 
by  a  twin  boundary  whose  shape  deformations  differ  by  a  simple  shear.  The 
original  lattice  is  termed  as  the  ‘parent’,  while  the  sheared  lattice  is  termed  as  the 
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‘twin’.  The  atomic  positions,  and  hence  the  corresponding  lattice  vectors 
between  these  atoms,  within  each  region  differ  by  a  finite  rotation,  typically 
either  a  reflection  or  180°  rotation  (Christian  &  Mahajan  1995),  although  more 
general  relationships  are  possible.  Nucleation  and  propagation  of  twins  are 
thought  to  take  place  by  one  or  more  mechanisms,  often  involving  the  formation 
and  motion  of  partial  dislocations  and  atomic  shuffles  needed  to  maintain  the 
orientation  relationships  between  the  twin  and  parent  (Bilby  &  Crocker  1965). 
Deformation  twinning  is  distinguished  from  plastic  slip  in  that  the  former  occurs 
by  collective  motion  of  defects,  resulting  in  a  quantized  amount  of  shear  that 
preserves  the  particular  orientation  relationship  between  the  twin  and  parent. 
Another  difference  is  that  twinning  is  polar  (i.e.  unidirectional)  while  usually  slip 
is  not.  Lattice  geometry  precludes  twinning  shears  of  equal  magnitude  and 
opposite  directions  on  the  same  plane,  while  typically  slip  may  occur  in  opposite 
directions  on  the  same  plane.  The  mechanical  work  of  deformation  twinning  is 
dissipative,  resulting  from  the  defect  motion  associated  with  shearing.  Possible 
energy  storage  is  associated  only  with  defects  left  behind  in  the  lattice,  for 
example  those  comprising  the  twin  boundary.  From  continuum  thermomecha¬ 
nics  considerations,  the  driving  force  for  twin  propagation  is  the  resolved  shear 
stress  in  the  direction  of  twinning  shear.  Detwinning,  i.e.  restoration  of  the 
twinned  lattice  to  its  original  orientation,  is  possible,  although  is  more  applicable 
to  phase  transformations  (Bhattacharya  1991)  and  less  applicable  to  mechanical 
twinning  of  the  kind  considered  here. 

Twins  are  usually  classified  as  type  I,  type  II  or  compound.  In  centro- 
symmetric  crystals,  the  lattice  vectors  in  the  twin  and  parent  for  a  type  I  twin 
are  related  by  either  a  reflection  in  the  habit  plane  or  rotation  of  180°  about  the 
direction  normal  to  this  plane.  For  a  type  II  twin,  the  lattice  vectors  are  related 
by  either  a  rotation  of  180°  about  the  shear  direction  or  a  reflection  in  the  plane 
normal  to  the  shear  direction.  The  orthogonal  transformation 


Q  = 


2m0®m0-l  (type  I), 
2s0<2>s0  —  1  (type  II) 


(2.1) 


relates  a  lattice  vector  a0  in  the  parent  to  a  lattice  vector  a  in  the  twin  via 
a  —  Qa0.  The  unit  normal  to  the  habit  plane  is  m0  and  the  direction  of  shear  is 
s0.  The  rank-two  unit  tensor  is  denoted  by  1.  For  a  compound  twin,  both 
definitions  in  (2.1)  are  equivalent. 


3.  Continuum  theory:  nonlinear  elasticity,  slip  and  twinning 

A  constitutive  framework  for  crystals  undergoing  large  thermoelastic,  plastic  and 
twinning  deformations  is  developed.  The  model  is  based  upon  principles  of 
continuum  mechanics  and  thermodynamics  of  single  crystal  behaviour,  as 
described  for  example  in  Clayton  (2005).  Plastic  deformation  is  limited  to  that 
resulting  from  dislocation  glide,  and  twinning  deformation  is  limited  to 
mechanical  twinning. 
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(a)  Kinematics 

A  local  volume  element  of  crystalline  material  is  identified  by  its  reference 
coordinates  X.  Let  x=<p(X,t)  denote  spatial  coordinates  of  the  element,  with  cp 
the  motion.  The  deformation  gradient  F  for  the  element  is  the  tangent  mapping 

F  =  T<px  =  dx/dX ,  (3.1) 

decomposed  multiplicatively  into  three  terms 

F  =  FeFIFp  =  FeF,  (3.2) 

where  FE  accounts  for  recoverable  thermoelastic  deformation  and  rigid  body 
rotation;  F1  =  J  1/3_F/  accounts  for  defect  kinematics  that  alter  the  lattice;  and 
Fp  accounts  for  lattice-preserving  plastic  slip.  Twinning  is  modelled  via  the 
isochoric  term  F1,  and  residual  volume  changes  associated  with  defects  are 
captured  by  the  Jacobian  determinant  J.  The  total  irreversible  deformation  is 
F  =  F 1  Fp . 

Introduced  next  are  sets  of  contravariant  and  covariant  vectors  denoting 
directions  and  planes,  respectively,  for  slip  and  twinning.  When  referred  to  the 
reference  lattice  prior  to  any  reorientation  by  twinning,  these  are  denoted  by 
{sq,  mj}  for  each  slip  system  i  and  {sJ0,mJ0}  for  each  twin  system  j.  The  total 
number  of  slip  and  twin  systems  is  denoted  by  n  and  w,  respectively.  Reference 
shearing  directions  and  plane  normals  are  all  of  unit  length,  and  each  pair  of 
contravariant  shear  direction  and  covariant  plane  normal  is  orthogonal 

so ' mo  =  0,  |«ol  =  \mb\  =  1  (V  *  = 

so '  mo  =  °,  lsol  =  lTOol  =  1  (Vj  =  1,  ••.,  w). 

During  the  course  of  twinning  deformation,  one  or  more  parts  (i.e.  twins)  of  the 
volume  element  of  crystal  undergo  a  rotation  relative  to  the  parent.  In  a  volume 
fraction  of  the  crystal  undergoing  twinning  via  mode  j,  the  slip  directions  and  slip 
plane  normals  transform  according  to  the  usual  rules  for  contravariant  and 
covariant  vectors,  i.e. 

sl)j  =  Q3s o,  m‘0j  =  ml0QjI 

(V  i  =  1, ...,  n  slip  systems;  V  j  =  1, ...,  w  twin  volumes),  (3.4) 

where  QJ  is  the  rotation  found  from  (2.1)  corresponding  to  particular  twin 
system  j.  For  example,  if  j  is  a  type  I  twin,  Q 3  —  2mJ0®m(  —  1,  while  if  j  is  a  type 

II  twin,  Q  =  2sq<8>Sq  —  1.  Note  from  (3.4)  that,  within  each  twinned  volume,  the 

updated  slip  directions  s  ^  and  slip  plane  normals  remain  orthogonal  and  of 
unit  length  for  each  i.  For  simplicity,  successive  twinning  is  not  considered, 
which  means  that  secondary  twins  do  not  form  within  already  twinned  regions 
and  reorientation  of  the  twinning  systems  {s30,  m]0}  is  not  considered.  The 
rotation  in  (3.4)  does  not  apply  to  the  volume  fraction  of  the  grain  comprising 
the  parent.  The  thermoelastic  deformation  affects  the  lattice  in  the  usual  manner 
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Figure  1.  Deformations  and  slip  system  geometry  for  crystal  deforming  by  elasticity,  slip 

and  twinning. 


of  crystal  plasticity  theory  (Asaro  1983) 

sl  =  FEsl0,  m1  =  (V  i  =  1, n  e  parent), 

s'j  =  FESQj,  m'j  =  m^F E~x  (V  *  =  1,  n  e  twins  j  =  1, w ),  (3.5) 

sJ  =  Fe s^,  m]  =  m\)FE~l  (V  j  =  1, w  G  parent). 

Figure  1  depicts  (3.1)-(3.5)  for  a  crystal  with  a  single  slip  system  and  a  single 
twin.  The  multiplicative  decomposition  (3.2)  implies  a  series  of  configurations 
of  the  material  element.  The  reference  configuration  is  labelled  as  B0  with 
coordinates  X ,  the  spatial  configuration  is  labelled  as  B  with  coordinates  x  and 
the  elastically  unloaded  intermediate  configuration  is  labelled  as  B .  Since  FE~ 1 
and  F  are  in  general  not  integrable,  continuous  coordinates  spanning  B  do  not 
exist  (Clayton  et  al.  2005).  However,  elastic  and  inelastic  deformations  act  as 
tangent  maps  between  configurations  via  FE  :  TB  — >  TB  and  F  :  TB0  —>  TB . 
In  this  illustration,  the  mappings  F  1  and  Fp  comprising  F  cannot  be  separately 
resolved  owing  to  the  effect  of  twinning  on  the  orientation  of  the  reference  lattice. 
The  spatial  velocity  gradient  following  from  (3.1)  and  (3.2)  is 

L  =  F  F~l  =FEFe~ 1  +  FeFIFi~1Fe~1 

+  FeFiFPFp~1Fi~1Fe~ 1  +  (1/3)  jj_1l.  (3.6) 


The  inelastic  velocity  gradient  referred  to  configuration  B  is 

F  F  =  L1  +  LP  +  (1/3)  Jj-1l, 

where 

W 

l1  ^V-^^/VVoK 

3= 1 

results  from  twinning  shears,  and 

n  w  /  n 

lp  =  (i-X-)X]^so0mo  +  Y1  [■fJ^Z^3sk'i^rnh 

i= 1  j= 1  \  i—l 


(3.7) 

(3.8) 


(3.9) 
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results  from  slip  in  both  the  parent  and  twinned  regions  of  the  crystal.  In  (3.8),  yJ 
is  the  predefined  shear  deformation  associated  with  twin  system  j.  In  (3.9),  y*  is 
the  slip  rate  on  system  i  in  the  parent  and  y)  is  the  slip  rate  on  system  i  within 
reoriented  twin  fraction  j.  The  volume  fraction  of  crystal  occupied  by  twin  j, 
measured  per  unit  volume  in  configuration  B ,  is  labelled  by  the  scalar  fJ>  0,  with 
rate  f  (Chin  et  al.  1969;  Van  Houtte  1978).  The  total  volume  fraction  of  twinned 
crystal  is  fx=sZ!fJ,  subject  to  the  restriction  0</r<l.  Detwinning  is  not 
considered;  hence,  f  J  >  0.  In  the  interior  summation  within  the  rightmost  term 
of  (3.9),  the  slip  directors  and  slip  plane  normals  in  the  twinned  regions  are  found 
according  to  (3.4),  where  the  particular  form  of  Q  corresponds  to  the  twin  with 
an  associated  value  of  / 3  in  the  outer  summation.  Note  that  since  for  each  slip  or 
twin  system,  the  shear  directions  and  plane  normals  are  orthogonal, 

jhtL1  =  0  =>  /  =  1,  /  JPtr {L1  +  Lp)  =  0  =>  JP  =  1,  (3.10) 

upon  assuming  that  at  t=  0,  F1—  Fp=  1  and  where  trA  —  Aaa  for  rank-two  matrix 
A.  Thus,  (3.8)-(3.10)  reflect  the  isochoric  nature  of  slip  and  twinning,  and 
together  with  (3.2)  and  (3.6)  require  that  all  volume  changes  be  accommodated 
thermoelastically  via  JE  and/or  by  defect  generation  via  J  ,  such  that 

J  =  y/det  gj det  G  det  F  =  JE  J  ,  (3-11) 

where  g  and  G  are  metric  tensors  associated  with  coordinate  systems  in 
configurations  B  and  B0,  with  respective  components  gab  =  dax-dbx  —  ga-  gb  and 
Gab  =  dAX-dBX  =  G\ ■  GB-  The  residual  volume  change  J  arising  from 
distributed  defects  such  as  dislocation  lines  follows  from  the  notion  that  in 
nonlinear  elasticity,  the  average  strain  of  a  body  containing  residual  stress  fields 
arising  from  internal  displacement  discontinuities  need  not  vanish  even  if  the 
traction  on  its  external  surfaces  vanishes  (Toupin  <fc  Rivlin  1960;  Wright  1982). 
Here,  only  volume  changes  are  considered,  e.g.  corresponding  to  random  defect 
distributions  imparting  no  preferred  directions  in  residual  strains,  although  more 
general  treatments  allowing  for  shape  changes  arising  from  the  local  elastic  stress 
fields  of  crystal  defects  have  been  posited  (Clayton  et  al.  2005).  In  §4e,  a 
relationship  between  J  and  the  dislocation  density  is  given. 

Two  non-dimensional  internal  state  variables  are  introduced  to  represent 
energetic  changes  associated  with  lattice  defects.  The  first  is  a  measure  of  the 
total  line  density  of  dislocations  that  accumulate  with  slip,  £  =  b^p^,  where  b  is  a 
scalar  Burgers  vector — or  a  constant  in  the  order  of  a  lattice  parameter  for 
crystals  containing  dislocation  families  with  different  Burgers  vectors — and  px  is 
the  total  length  of  such  dislocations  per  unit  volume  in  B.  The  second  is  a 
measure  of  the  total  area  density  of  twin  boundaries,  £  =  \/brj T,  where  is  the 
total  area  of  twin  boundaries  per  unit  volume  in  B . 

( b )  Constitutive  assumptions 

Let  p  and  p0  denote  the  mass  density  of  the  solid  in  current  and  reference 
configurations,  respectively,  related  by  p0= pj.  Let  p—  pJE  —  p0J  denote  the 
mass  density  in  configuration  B .  The  Helmholtz  free  energy  per  unit  volume  in  B 
is  W  =  p\p,  where  is  the  free  energy  per  unit  mass.  The  free  energy  exhibits 
the  dependencies 

w  =  w{EE,e^xAfJ}),  (3.12) 
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where  6  is  the  absolute  temperature;  %  and  £  are  respectively,  related  to  the 
densities  of  dislocations  and  twin  boundaries  as  discussed  in  §3 a;  the  set  {/J} 
includes  each  of  the  twin  fractions,  and 

=\{C%-  5ap)  =  i  (. FEaagabF%b  -  5af})  (3.13) 

is  a  finite  elastic  strain  tensor  associated  with  the  covariant  elastic  deformation 
tensor  CE,  with  Cartesian  metric  8ap  used  in  configuration  B . 


(c)  Thermodynamics 

Local  forms  of  the  balance  of  energy  and  dissipation  inequality  referred  to  the 
reference  configuration  are  written  respectively,  as 

p0e  =  2  :  E  -  V0- Q  +  p0r,  2  :  E  -  p0(\p  +  pd)  -  ^ Vo0- Q>  0.  (3.14) 

Here,  e=\p  +  6p  is  the  internal  energy  per  unit  mass,  with  p  the  entropy  per  unit 
mass.  The  symmetric  second  Piola-Kirchhoff  stress  2  is  related  to  the  first  Piola- 
Kirchhoff  stress  P  and  the  Cauchy  stress  a  by 

77?  —  1A  t)ClB  j 77?  —  IA  ^.ab  77? — 1 B  /o  1  r\ 

E  =  F.a  P  =  JF.a  a  F.b  .  (3.15) 

Also,  E={1/2)FTF-  G  is  the  right  Cauchy-Green  strain;  V0  denotes  a 
covariant  derivative  on  B0 ;  Q  is  the  heat  flux  vector;  and  r  denotes  other  heat 
sources.  The  stress  power  per  unit  intermediate  volume  can  be  written  as 


J  l2ABEAB  =  (  JEFB~W"ab' 


)  {FE~laLedFf)  =  MfL%  (3.16) 


J  _au  77? T 

&  Qca^.c 

where  M=  CE2,  2  =  JE FE~loFE~T ,  and 

L  =  Fe~1LFe  =  Fe~1FE  +  L1  +  LP  +  (1/3)  JJ 
The  rate  of  free  energy  per  unit  intermediate  volume  is 


1. 


v= am  =  J 


7  -1 


-1A 


Poi^-E.AF.a  '  \j/ 


and  following  from  (3.17)  and  the  symmetry  properties  of  2  and  EJ 


M TE% 


2>uEsy  +  M'a  Bp  +M'a  LPp  +  (1/3)  J  J  ~lM'p 
Using  the  chain  rule  to  expand  the  rate  of  W  of  (3.12), 

W  = 


dW  ■  E  dW  ■  dW  ■  dW  ■  dW -j 

dE^:E  +^6d  +  ^+^  +  Wf  ’ 


(3.17) 

(3.18) 

(3.19) 

(3.20) 


with  summation  implied  over  w  twin  fractions  j,  the  entropy  inequality  in  (3.14) 
can  be  written  as 


f -  dW\ - 
dW  ■_  dW  ■_  dW 


:(i/  +  iP  +  3j1) 
fj-^V0-q>  0, 


(3.21) 
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where  N  =  pq  is  the  entropy  per  unit  intermediate  volume;  V6=V06F  1  is 
the  intermediate  temperature  gradient;  q  =  J  F  Q  is  the  intermediate  heat 
flux;  and 

n  =  M-m  (3.22) 

is  a  version  of  Eshelby’s  energy-momentum  tensor  (Eshelby  1975)  mapped  to 
configuration  B .  Following  standard  arguments  (Coleman  &  Noll  1963),  stress- 
elastic  strain  and  entropy-temperature  relationships  are  deduced  from  (3.21)  as 


-  3W 
^~~3Ee' 


3W 


3W 


N=-^~,  a  =  Je~1Fe - = 

36  ’  3Ee 


-iET 


V  =~P 


--1 


3W 
~36  ' 


(3.23) 


The  final  term  on  the  left-hand  side  of  (3.21)  can  be  made  non-negative  upon 
prescription  of  the  conduction  law 

q=-KW,  -V6-q=V6-KW>  0,  (3.24) 

where  K  is  a  symmetric  and  positive  definite  matrix  of  thermal  conductivity. 
Applying  (3.23)  and  (3.24), 


77 


L1  +77 


^  +  +  <-> 

is  the  reduced  dissipation  inequality.  In  the  absence  of  temperature  gradients, 
(3.25)  requires  that  the  energy  dissipated  by  twinning,  slip  and  residual  volume 
changes  must  exceed  the  rate  of  energy  storage  associated  with  defects, 
specifically  dislocations  and  twin  boundaries.  From  (3.8)-(3.10),  the  energies 
dissipated  from  twinning  and  slip,  respectively,  can  be  written  as 


77 


L1  = 


J2  V,  n  lp  =  ( 1  -  /T)  J]  f'y  +  XI  Y, 


(3.26) 


3=1 


i  =1 


3=1 


i  =1 


where  the  driving  forces  are  resolved  stresses  on  each  slip  or  habit  plane,  acting 
in  the  direction  of  shear 


?  ja  -b  7 

y  mJb, 


=  JV 


G a 


=  jeP 


a  n  m. 


W- 


(3.27) 


The  specific  heat  capacity  at  constant  elastic  strain,  measured  per  unit  volume  in 
B ,  is  introduced  as 

c  =  3E  jdd  =  -6d2W/d6\  (3.28) 

where  (3.23)  has  been  used  and  E  =  pe  is  the  internal  energy  per  intermediate 
volume.  Multiplying  the  first  of  (3.14)  by  J  1 ,  and  using  (3.15)-(3.20),  (3.22)-(3.24) 
and  (3.28),  the  energy  balance  can  be  written  as 


cd  =  77  :  (7/  +  LP )  + 


^tr77  +  6 


J 
3 J 


-6/3  :  EE 


+  V-KV8 


_'3W  _  d2W\ 

“  W~ffd6dfJ 


3W 


d2W\, 

36dZjl 


3W 


3P 


— — V 

36 dp  J 


(3.29) 


where  non-mechanical  sources  rin  (3.14)  are  assumed  absent  in  (3.29)  and  hereafter. 
The  stress-temperature  coefficients  are 


0  =  - 32W/363Ee , 


(3.30) 
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and  the  notation  for  the  heat  conduction  term  in  (3.29)  is  Va  =  Va  +  J  1  ( J  F  a1A)  A, 
with  Va  =  VoaF  a  .  Notethat  Va  =  Va  when  the  compatibility  conditions  F  a[A,B]  =  0 
apply,  in  which  case  ( dJ  /dFaA)  A  =  0. 


( d )  Representative  free  energy  potential 


A  more  specific  form  of  (3.12)  is  now  specified  for  anisotropic  crystals  that  may 
undergo  large  elastic  deformations.  The  free  energy  is  decomposed  additively  as 

w  =  WE(Ee,  d,  {f})  +  Y(0)  +  WR(^  C,  0),  (3.31) 

where  W  accounts  for  the  thermoelastic  response;  Y  accounts  for  the  purely 
thermal  energy;  and  W  accounts  for  residual  energy  of  lattice  defects.  The 
thermoelastic  energy  consists  of  three  terms 

WE  0Q),  (3.32) 

with  the  first  term  of  (3.32)  accounting  for  materially  linear,  but  geometrically 
nonlinear,  mechanical  effects;  the  second  accounting  for  nonlinear  elastic  effects 
important  at  high  pressures  (Brugger  1964;  Graham  &;  Brooks  1971);  and  the 
third  accounting  for  thermoelastic  coupling.  Here,  80  is  a  constant  temperature  at 
which  the  lattice  parameters  exhibit  their  reference  lengths,  and  the  remaining 
coefficients  in  (3.32)  consist  of  partial  derivatives  of  the  free  energy  at  null 
elastic  strain 


CaW  = 


d2WE 


dEEpdFEs 


C 


afixSerl)  _ 


3)7 rE 


dAW 


Ee=  0 


dE&E&E% 


,  = 


ee= o 


d2WE 


Ee= 0 


(3.33) 

Coefficients  in  (3.33)  may  depend  on  temperature;  when  measured  at  a  particular 
temperature,  these  are  referred  to  as  isothermal  elastic  constants.  In  anisotropic 
solids,  the  coefficients  in  (3.33)  depend  upon  the  orientation  of  the  lattice  in 
configuration  B .  When  twinning  is  involved,  the  orientations  of  the  original 
reference  lattice  (parent)  and  each  twinned  region  differ.  It  is  assumed  that  the 
elastic  deformation  FE  and  strain  Ee  act  uniformly  over  the  parent  and  twins 
comprising  each  volume  element.  The  energy  (3.32)  is  thus  partitioned  into 
contributions  from  the  parent  and  each  twin 

VE  =  \  (1  -  fT)EEyCfxSEEs  +  i  (1  - 

-  (1  -  fr)P ?E%(6  -  d0)  +  J2  f  (2  E«^tx6EEs 

+  \E%C^E*sE%-p%E*(t{e-e^  (3.34) 

where  Cg^*5,  Q“^x<5f<^  and  /3g^  refer  to  coefficients  for  the  parent  lattice,  and  where 
for  each  twin  j  (Kalidindi  1998), 

/p“/?x<5  _  nJan  JP  n  Jx  n  ^Px^  _  r>Pa  n  tP  n  tx  /g  je  ^  j<j> 

^ jT  —  .£  V  4  H  .ip  h?  .yM)  )  —  W  .ipW  .yW  .11  W  .l  W  .kW  .X  '*-'0  ) 

(S?t  =  Q’"xQ  (3.35) 
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The  stress-strain-temperature  relationships  followed  from  (3.23)  and  (3.31)— 
(3.34)  are 

2#  =  C apxSEE6  +  Ic Eex5E% -Pa\e-  do),  (3.36) 

where  the  effective  coefficients 

W  W 

Cafix5  =  (l  -  /T)  cfxi  +  J2  f£jTX5,  Caft^  =  (1  -  /t)Co  W  +  J2  f^TX5^, 

j= 1  '  j= 1 

W 

Paft  =  (1  -  fi)P?  +  (3-37) 

l=i 

From  (3.37),  the  rationale  for  inclusion  of  twin  fractions  in  free  energy  functions 
(3.12)  and  (3.31)  is  now  clear:  the  effective  thermoelastic  moduli  depend  on 
the  evolving  twin  fractions.  The  associated  rate  of  thermoelastic  free  energy 
change  is 


dWE 

w 


/■'  =AJ 


A,  =\ee 


C“ 


a(3\  ^jT 


w  771-& 

^  ^  ^  £(P 


ofrc5  _  pB 


0 

0L& 


) 


+  5* 

E 


/p1 

apV^jT 


aPxSecj) _ 0a/?x<5#  ) 


0 


) 

(3.38) 


with  summation  applied  over  j.  Consider  a  situation  in  which  strains  EEq  = 
aap(d  —  6q)  arise  from  temperature  change.  The  following  relationship  emerges 
between  the  thermal  stress,  thermal  expansion  (ciap)  and  elasticity  coefficients 

r#  =  0“”%  +  l-c^“*axld,t(e- «„).  (3.39) 

Note  that  Ca/3*5  =  C(“Wz5),  £****  =  and  j}<4  =  qW)  foUow  auto_ 

matically  from  (3.33).  The  Voigt  notation  (Brugger  1964)  exploits  the  symmetry 
of  the  elastic  stress  2,  elastic  strain  Ee  and  these  coefficients,  with  pairs  of  indices 
11  —>  1, 22— >2, 33— >3, 23— >4, 13— >5 and  12— >6 and 2af}  ->  2A,  2EEfi  ->  Eea (1  +  <5^), 
with  Greek  indices  (a,  0=1,2,  3)  and  capital  indices  {A  =  \,  2,  ...,  6).  Relation  (3.36) 
then  can  be  expressed  compactly  as 


SA 


Cab^e  + 


1 

2 


CabcEbEe  ~ @  ~  0o), 


(3.40) 


where  Cab  and  CAbc  are  the  second-  and  third-order  elastic  coefficients  and 
summation  proceeds  over  duplicate  covariant  indices.  The  thermal  energy  in  (3.31)  is 
prescribed  as  (Clayton  2005) 

Y  =  -cd\n(d/e0).  (3.41) 

Finally,  the  residual  energy  of  (3.31)  is  specified  as 

wr = ^(Kle + *2c2 + K3e?),  (3.42) 


where  p  is  an  elastic  shear  modulus  that  may  depend  on  temperature  and  /q,  k2  and  /c3 
are  dimensionless  constants  that  scale  the  energies  associated  with  each  internal 
variable.  Recalling  from  §3a  that  £  =  the  first  term  on  the  right-hand  side  of 
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(3.42)  provides  a  linear  dependence  of  residual  energy  on  the  line  density  of 
dislocations,  following  Clayton  (2005)  and  references  therein.  Simple  arguments  then 
show  that  2W^  =  Ki)xb2 ,  where  W^>  is  the  total  energy  per  unit  line  length  of 
dislocations,  including  self-  and  interaction  energies,  core  energy  and  stacking  fault 
energy  if  the  dislocations  are  partial.  Recalling  that  £  =  y/brj T,  the  second  term 
provides  a  linear  dependence  of  residual  energy  on  the  area  per  unit  volume  of  twin 
boundaries  ryx-  Similarly,  2  Wr  =  k^/a  b,  where  Wr  is  the  twin  boundary  energy  per  unit 
area.  Often,  2  ITt~  WSF,  where  W§f  is  the  intrinsic  or  extrinsic  stacking  fault  energy 
(Hirth  &  Lothe  1982;  Bernstein  &  Tadmor  2004).  The  third  term  accounts  for 
interaction  energies  between  twin  boundaries  and  dislocations.  Alternative  forms  of 
(3.42)  are  of  course  possible  and  can  be  introduced  as  needed  for  particular  materials. 
The  rate  of  temperature  increase  (3.29)  can  now  be  written  as 

ce  =  p'Wp  +  e{d^^L-p-.  ee\  +  v-kvo,  (3.43) 


where 


P'  =  I' 


fi  —  6 


d/i 

Jo 


Vt 

Vpt 


Pt  + 


*,63/2 


Pt 

VJt. 


Vt 


+  ( A,-ed^j)fJ 


3 


do 


w 


P-1 


(3.44) 


is  often  called  the  Taylor-Quinney  parameter  (Taylor  &  Quinney  1934),  such  that 
1  —  J  is  the  ratio  of  inelastic  stress  power  Wp  =  77  :  (L1  +  LP  +  (1/3)  JJ  xl) 
converted  to  residual  elastic  energy  in  the  lattice. 


(e)  Kinetics 

Here,  rate-dependent,  i.e.  viscoplastic,  relationships  describe  shearing  associ¬ 
ated  with  slip  (Hutchinson  1976;  Asaro  1983)  and  twinning  (Kalidindi  1998). 
These  kinetic  equations  read  as 


y)  =  is 


7  T 

(f3) 

j3 

J 

(3.45) 


In  (3.45),  an(i  7 t  are  material  parameters  with  dimensions  of  1  /£;  to  and  p  are 
dimensionless  parameters;  and  <f  and  g3  are  evolving  resistances — positive 
scalars  with  dimensions  of  stress — to  deformation  in  the  parent  by  slip  on  system 
i  and  twinning  on  system  j,  respectively.  In  the  second  of  (3.45),  g*  denote 
resistances  on  slip  system  i  within  reoriented  twin  j.  In  the  third  of  (3.45), 
2 (rJ)  =  f3  +  \t3\.  From  (3.26)  and  (3.45),  these  rates  are  always  dissipative  since 
for  each  slip  system,  r*y*>0  and  r)y)  >  0  and  for  each  twin,  /J(fJ)yJ  >  0. 
Because  f  J  =  0  for  f3  <  0,  the  unidirectional  nature  of  twinning  is  respected.  In 
the  limit  m— ►  °°  or  p— ►  oo5  rate-independent  behaviour  is  attained.  In  a  general 
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sense,  slip  and  twin  resistances  evolve  as 

$  =  9(EE,  o,  ?,  C,  {f},  {g}),  -g)  =  %(EE,  9,  C,  {f},  {$}), 

^  =  ^(EE,0,^,{f},{g}),  (3.46) 

where  the  rates  depend  not  only  on  the  set  of  state  variables  that  explicitly  enter 
the  free  energy  (3.12),  but  also  on  the  set  of  hidden  variables  {g}=  [g\  g3,  Jfj} 
with  i=l,...,n  and  j=l,  ...,  w.  Evolution  equations  for  f  and  £,  reflecting 
densities  of  dislocations  and  twin  boundaries,  complete  the  model.  Analogous  to 
hardening  in  (3.46),  generic  evolution  equations  are  written  as 

|  =  t(EE,  9 ,  ?,  £,  {/},  {g}),  t  =  t(EE ,  9,  ?,  ?,  {f},  {g}).  (3.47) 

Possible  impedances  of  slip  or  twinning  via  slip-slip  interactions,  slip-twin 
interactions  and  twin-twin  interactions  entering  (3.46)  may  depend  in  a  complex 
manner  on  a  number  of  factors,  including  geometric  relationships  between 
interacting  systems,  temperature,  crystal  structure  and  defect  content  (Christian  &; 
Mahajan  1995;  Kalidindi  1998;  Castaing  et  al.  2002;  Wu  et  al.  2007).  The  theory  is 
applicable  over  a  potentially  wide  temperature  range  via  dependence  of  initial 
values  of  resistances  to  slip  and  twinning  (i.e.  initial  values  of  {<)}  in  kinetic 
equations  (3.45))  on  temperature;  in  addition,  these  variables  evolve  with  the 
history  of  deformation  and  temperature  according  to  (3.46). 


4.  Application:  sapphire  single  crystals 

The  framework  of  §3  is  specialized  to  describe  sapphire.  Crystal  structure, 
physical  properties,  initial  yield  mechanisms,  strain  hardening,  defect  accumu¬ 
lation  and  residual  volume  changes  are  considered.  Whenever  possible, 
experimental  data  motivate  aspects  of  the  model  and  provide  associated  parameters. 

(a)  Crystal  structure  and  properties 

The  atomic  structure  of  sapphire  is  depicted  in  terms  of  the  Al  nuclei  in 
figure  2,  following  Kronberg  (1957).  The  fundamental  unit  cell  is  rhombohedral 
(lattice  parameter  ao  =  0.512  nm,  bond  angle  55.3°).  The  hexagonal  unit  cell, 
while  consisting  of  more  atoms,  is  convenient  for  describing  mechanical 
behaviour.  Two  types  of  hexagonal  cells  are  encountered  in  the  literature.  One 
is  the  morphological  cell  shown  in  figure  2  a,  consisting  of  three  layers  of  Al 
cations  each  two-third  full,  with  a  close-packed  full  layer  of  O  anions  in 
between  each  layer  of  Al  cations.  The  Al  nuclei  occupy  octahedral  interstitial 
sites  between  each  hexagonal  layer  of  O  atoms  (Kronberg  1957).  The 
morphological  hexagonal  lattice  parameters  are  A  =  0.475nm  and  Cm  = 
0.649  nm.  Figure  2b  is  of  the  structural  cell,  consisting  of  a  second  stack  of 
each  of  the  three  layers  of  Al  cations  and  O  anions,  differing  from  the  morpho¬ 
logical  unit  cell  by  a  rotation  of  180°  about  the  c-axis  [0001].  The  lattice 
parameters  of  the  structural  hexagonal  cell  are  A/  =  0.475nm  and  Cs  =  2Cm  = 
1.297  nm,  as  shown  in  figure  2.  Elastic  constants  are  listed  in  table  1  (in  the 
electronic  supplementary  material);  because  sapphire  has  trigonal  symmetry, 
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rotations  of  180°  about  the  c-axis  do  not  fall  within  the  point  group  of 
symmetry  operations  for  the  crystal.  Hence,  the  signs  of  C14  and  several  of  the 
third-order  constants  depend  upon  the  choice  of  unit  cell  (Winey  et  al.  2001). 
In  table  1  and  henceforth,  all  properties  are  referred  to  the  structural  unit  cell. 
The  full  matrices  of  anisotropic  elastic  coefficients  in  (3.40)  are  populated  from 
the  six  independent  second-order  coefficients  and  14  independent  third-order 
coefficients  using  the  scheme  of  Teodosiu  (1982),  for  example.  Relevant  bulk 
physical  properties  with  supporting  references  are  listed  in  table  2  (in  the 
electronic  supplementary  material). 

Prominent  slip  and  twin  systems  are  listed  in  tables  3  and  4  (in  the  electronic 
supplementary  material),  respectively,  and  are  illustrated  in  figure  2.  While  a 
number  of  different  families  of  pyramidal  systems  are  crystallographically 
conceivable  (Snow  &;  Heuer  1973),  quantitative  yield  data  seem  scarce  except  for 
the  <1101  >  {1011}  systems,  with  available  data  limited  to  tensile  deformation  of 
whiskers  at  temperatures  in  excess  of  2000  K  (Tressler  &;  Barber  1974).  While 
basal  and  prism  slip  can  occur  in  either  direction  for  a  given  system,  pyramidal 
slip  is  thought  unidirectional  (Heuer  et  al.  1998)  and  only  occurs  when  the 
resolved  shear  stress  acts  in  a  positive  sense  with  respect  to  the  c-axis,  e.g.  tensile 
loading  along  [0001].  Twinning  is  also  unidirectional;  for  example,  rhombohedral 
twinning  occurs  only  when  the  resolved  shear  stress  acts  in  a  negative  sense  with 
respect  to  the  c-axis,  e.g.  compressive  loading  along  [0001].  Thus,  pure  tensile 
deformation  along  [0001]  can  only  be  accommodated  by  pyramidal  slip,  elasticity 
or  fracture,  while  pure  compression  along  [0001]  can  only  be  accommodated  by 
rhombohedral  twinning,  elasticity  or  fracture.  The  slip  directions  for  uni¬ 
directional  mechanisms  are  tabulated  here  such  that  the  resolved  shear  stress 
for  a  given  slip  or  twin  system  must  be  positive  to  enact  shear  on  that  system. 
Atomistic  simulations  of  hypervelocity  impact  of  sapphire  have  predicted  basal 
and  pyramidal  slip  and  basal  and  rhombohedral  twinning  (Zhang  et  al.  2007). 
Bourne  et  al.  (2007)  observed  prism  and  basal  dislocations,  twins,  cleavage 
fracture  and  grain  boundary  fracture  in  specimens  recovered  from  impact 
experiments  on  polycrystalline  alumina,  with  activity  or  inactivity  of  certain 
mechanisms  dependent  on  the  impact  stress. 

Quantification  of  resistances  for  slip,  twinning  and  fracture  in  sapphire 
remains  an  area  of  active  research  (Tymiak  &  Gerberich  2007).  Fundamental 
arguments  based  on  elasticity  theory  and  lattice  periodicity  provide 
order-of-magnitude  estimates  for  critical  shear  stresses  that  will  be  refined 
upon  examination  of  experimental  data.  The  isotropic  Peierls  stress  (Peierls 
1940)  has  been  deemed  useful  for  describing  basal  and  prism  slip  in  sapphire 
(Farber  et  al.  1993) 


TC 


K  = 


1—v  (edge  dislocation) 
1  (screw  dislocation) 


(slip), 


(4.1) 


where  d  is  the  spacing  between  slip  planes  and  v  is  Poisson’s  ratio.  Eshelby  (1949) 
derived  expressions  for  anisotropic  crystals;  however,  owing  to  their  complexity, 
the  need  for  consideration  of  the  orientation  of  the  dislocation  line  with  respect  to 
the  lattice  vectors  and  limitations  of  the  sinusoidal  interatomic  potential  used 
in  the  derivations,  the  anisotropic  solutions  are  usually  not  considered,  and  are 
not  pursued  here.  Table  5  (in  the  electronic  supplementary  material)  lists  the 
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Peierls  stresses  for  sapphire;  note  that  partial  Burgers  vectors  are  involved  for 
prism  slip,  leading  to  a  lower  Peierls  stress  and  a  preference  over  basal  slip  at  low 
temperatures  (Lagerlof  et  al.  1994),  although  the  original  treatment  of  Peierls 
(1940)  did  not  consider  partials.  For  these  computations,  the  rhombohedral  shear 
modulus  at  room  temperature  is  used  as  an  estimate  (table  2),  along  with  the 
effective  Poisson’s  ratio  for  alumina,  ^  =  0.24  (Holm  et  al.  1999).  The  expression 
below  has  been  postulated  for  twin  nucleation  (Hirth  &  Lothe  1982;  Lagerlof 
et  al.  1994) 

tc  =  WsF/b  (twin  nucleation) ,  (4.2) 

recalling  that  TFsp  is  the  stacking  fault  energy  associated  with  partial 
dislocations  of  Burgers  vector  b  responsible  for  twinning  on  the  system  of 
interest.  Table  6  (in  the  electronic  supplementary  material)  lists  corresponding 
values  for  nucleation  of  basal  and  rhombohedral  twins.  A  theoretical  first-order 
estimate  for  shear  fracture  is  (Frenkel  1926) 

tc  =  iab/(2nd)  (shear  fracture) ,  (4.3) 

with  values  listed  in  table  7  (in  the  electronic  supplementary  material)  for 
possible  fracture  planes  in  sapphire.  Experiments  indicate  that  the  preferred 
cleavage  planes  are  {1010}  prism  planes  and  {1012}  rhombohedral  planes 
(Wiederhorn  1969;  Azhdari  &;  Nemat-Nasser  1998).  The  prism  planes  listed  in 
table  7  are  of  the  type  {1010},  and  hence  comprise  a  different  family  than  the 
{1120}  prism  planes  considered  for  slip  in  tables  3  and  5  and  figure  2. 
Approximation  (4.3)  does  not  account  for  electrostatic  forces  in  ionic  crystals 
apart  from  their  influence  on  the  shear  modulus,  and  hence  neglects  the  strong 
Coulomb  forces  contributing  to  the  high  fracture  energy  of  the  {0001}  planes, 
relative  to  the  prism  and  rhombohedral  planes,  in  sapphire  (Wiederhorn  1969). 
Experimental  measurements  of  cohesive  energy  or  cleavage  strength  of  {10ll} 
pyramidal  planes  are  not  available.  Owing  to  sapphire’s  brittle  nature  at  low 
temperatures,  measurements  of  the  yield  mechanisms  corresponding  to  (4.1)- 
(4.3)  must  occur  at  high  pressures  that  suppress  tensile  fracture,  e.g.  indentation 
(Tymiak  &;  Gerberich  2007)  or  confined  compression  (Graham  Sz  Brooks  1971; 
Castaing  et  al.  1981;  Scott  &  Orr  1983;  Lankford  et  al.  1998).  Graham  Sz  Brooks 
(1971)  suggested  that  shear  failure  via  attainment  of  the  theoretical  strength 
(4.3)  could  occur  in  shock-loaded  sapphire. 

( b )  Thermo  elasticity  and  critical  shear  stresses  in  uniaxial  strain 

An  analysis  of  the  thermoelastic  response  of  sapphire  under  adiabatic  uniaxial 
strain  boundary  conditions  enables  consideration  of  the  behaviour  of  oriented 
crystals  up  to  the  yield  point  or  HEL.  Driving  forces  for  slip,  twinning  and 
fracture  can  be  computed  and  compared  with  (4.1)-(4.3)  and  observations  from 
shock  physics  experiments  (Barker  &  Hollenbach  1970;  Graham  Sz  Brooks  1971; 
Chen  Sz  Howitt  1998;  Reinhart  et  al.  2006)  to  provide  more  suitable  estimates  for 
yield  and  failure  criteria.  This  analysis  also  provides  insight  into  which  inelastic 
deformation  mechanisms  are  most  likely  responsible  for  macroscopic  yielding  in 
shock  physics  experiments.  In  conventional  plate  impact  experiments,  the 
sequence  of  activity  of  yield  and  fracture  mechanisms  can  usually  only 
be  inferred  indirectly  (Bourne  et  al.  2007),  since  visual  characterization  of 
defects  (e.g.  microscopy  of  dislocations,  twins  and  micro-cracks)  is  generally  not 
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Figure  3.  Axial  Cauchy  stress  for  single  crystals:  elasticity  predictions  and  shock  physics 
experiments  (solid  line,  0°  nonlinear  elasticity;  dotted  line,  60°  nonlinear  elasticity;  short-dashed 
line,  90°  nonlinear  elasticity;  dot-dashed  line,  rhombohedral  nonlinear  elasticity;  long-dashed 
line,  0°  linear  elasticity;  white  squares,  Barker  &  Hollenbach  (0°);  black  squares,  Graham  &  Brooks 
(0°,  60°  and  90°);  white  circles,  Reinhart  et  al.  (0°)). 

possible  in  situ ,  but  rather  is  conducted  on  recovered  specimens  that  fracture  on 
tensile  release,  if  not  during  the  impact  process.  Stress-induced  phase 
transformations  are  thought  to  occur  only  at  very  high  pressures  (approx. 
79  GPa;  Reinhart  et  al.  2006)  and  are  not  considered. 

For  uniaxial  straining  along  the  Cartesian  three-axis,  the  deformation  gradient 
(3.1)  is  F=  1  +  (A  V/  V0)g3  <£>  G3,  where  AV=  V—  V0  is  the  volume  change  and 
Vo  is  the  reference  volume  of  the  crystal.  Thermoelastic  responses  are  computed 
for  crystals  oriented  with  respect  to  the  loading  axis  at  0°  (along  [0001]),  60° 
(along  [1123])  and  90°  (along  [1210]),  and  with  the  loading  axis  normal  to  a 
rhombohedral  plane  (along  [1012]).  Shown  in  figure  3  is  the  axial  stress 
component  a  =  a'’"’  computed  for  each  orientation,  accounting  for  anisotropic 
nonlinear  elasticity  and  thermal  expansion.  Also  shown  is  the  response  for  a 
0°-oriented  crystal  computed  using  only  the  second-order  elastic  constants,  i.e. 
Cabc~ 0  in  (3.40),  and  labelled  ‘linear  elasticity’.  Experimental  data  points  are 
provided  from  the  aforementioned  sources,  and  dotted  horizontal  lines  indicate 
maximum  observed  HEL  values  quoted  by  Graham  &  Brooks  (1971).  Among  the 
nonlinear  elastic  predictions,  minute  differences  in  stress  arise  between  crystals  of 
different  orientations  below  the  HEL.  Among  corresponding  experiments, 
differences  among  crystals  of  different  orientations  are  not  discernable  in  the 
elastic  regime  (Graham  &;  Brooks  1971).  From  figure  3,  nonlinear  elasticity 
provides  an  adequate  match  to  the  experimental  results  for  V/  Vo  >0.95,  while 
linear  elasticity  is  inadequate  for  volume  changes  V/  Vo  <0.98.  For  each 
orientation,  temperature  rises  predicted  by  the  nonlinear  theory  from  an  initial 
temperature  of  0n = 300  K  were  of  the  order  of  20  K  at  V/  V0  =  0.95.  In  the 
experiments,  significant  stress  drops  occur  for  V/ Vo  <0.95;  comparison  with 
the  hydrostat  of  sapphire  indicates  a  significant  loss  of  shear  strength  above  the 
HEL,  pointing  to  possible  fracture  processes.  Values  of  the  post-HEL  shear 
strength  have  been  reported  in  the  range  of  4. 0-6.0  GPa  (Graham  &  Brooks 
1971;  Munson  &;  Lawrence  1979;  Reinhart  et  al.  2006). 
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Figure  4.  Ratio  of  maximum  shear  stress  computed  from  nonlinear  elasticity  to  theoretical 
estimates  of  Peierls  stress,  shear  fracture  strength  or  twin  nucleation  stress:  (a, b)  60°  orientation; 
dotted  lines,  HEL  (Graham  &  Brooks  1971)  ( c,d )  rhombohedral  orientation;  dotted  lines,  HEL 
(Hare  et  al.  2002).  ( a,c )  solid  line,  basal  edge  Peierls;  dotted  line,  basal  screw  Peierls;  dashed  line, 
prism  edge  Peierls;  dot-dashed  line,  prism  screw  Peierls;  (b,d)  solid  line,  basal  theoretical  strength; 
dotted  line,  prism  theoretical  strength;  small-dashed  line,  pyramidal  theoretical  strength;  double 
dot-dashed  line,  rhombohedral  theoretical  strength;  long  dashed  line,  basal  twinning;  single 
dot-dashed  line,  rhombohedral  twinning. 


For  each  orientation,  during  the  uniaxial  loading  process,  resolved  shear 
stresses  r  for  various  inelastic  deformation  modes  were  computed.  For  slip, 
r=max|f*|,  where  i€=all  basal  or  prism  slip  systems;  for  twinning,  t—  maxfJ, 
where  all  basal  or  rhombohedral  twin  systems.  For  fracture,  r  is  computed  as 
the  maximum  resolved  Kirchhoff  stress  in  any  in-plane  direction  among  each  of 
the  basal,  prism,  pyramidal  or  rhombohedral  planes.  Pyramidal  slip  is  not 
considered  in  this  analysis,  since  experiments  are  inconclusive  regarding  its 
occurrence  in  the  room  temperature  regime.  Pyramidal  slip  has  been  observed  in 
c-axis  sapphire  filaments  stretched  at  temperatures  in  excess  of  2000  K,  and  is 
thought  to  be  controlled  by  point  defect  diffusion  and  dislocation  climb  (Tressler  & 
Barber  1974).  Selected  results  are  shown  in  figure  4,  with  the  shear  stresses  r 
normalized  by  the  theoretical  estimates  tq  for  slip,  twinning  and  fracture  of  (4.1)- 
(4.3)  and  tables  5-7.  Volumetric  strains  corresponding  to  the  HEL  (Graham  & 
Brooks  1971;  Hare  et  al.  2002)  are  shown  with  vertical  dotted  lines.  In  the 
60°-oriented  crystal,  as  shown  in  figure  4a,  prism  slip  could  be  activated  well  before 
attainment  of  the  HEL,  if  the  Peierls  approximation  (4.1)  is  accurate.  This  agrees 
with  experimental  observations  of  l/3<ll00>  partial  dislocations  in  shock- 
loaded  crystals  of  the  same  orientation  (Chen  &  Howitt  1998).  Stresses  that  would 
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cause  twinning  or  shear  fracture  are  depicted  in  figure  4 h  for  a  60°-oriented  crystal. 
Among  the  mechanisms  shown  in  figure  4  b,  basal  twinning  is  predicted  as  most 
probable,  in  agreement  with  the  experiments  (Chen  Sz  Howitt  1998).  Regarding 
fracture,  failure  on  rhombohedral  or  prism  planes  is  most  likely,  but  theoretical 
fracture  stresses  from  (4.3)  are  not  attained.  Shown  in  figure  4c  are  shear  stresses 
for  glide  normalized  by  the  corresponding  Peierls  barriers  (4.1)  for  a  crystal  in  the 
rhombohedral  (R)  orientation.  Prism  slip  seems  feasible  here  at  or  below  the  HEL. 
Stresses  for  twinning  or  fracture  are  shown  in  figure  4  d.  Basal  twinning  appears 
most  likely  among  these  mechanisms,  and  normalized  driving  forces  for  shear 
fracture  on  pyramidal,  prism  and  basal  planes  are  also  significant.  Similar  analysis 
of  results  for  the  90°  orientation  (not  shown  in  figure  4)  predicts  that  for  a  90°- 
oriented  crystal,  prism  slip  below  the  HEL  is  possible  according  to  the  Peierls 
barrier,  twinning  is  unlikely  and  pyramidal  fracture  is  the  most  likely  shear  failure 
mechanism.  In  the  0°-oriented  crystal  (not  shown  in  figure  4),  slip  is  impossible  for 
the  present  set  of  loading  conditions.  Rhombohedral  twinning  is  most  likely, 
and  could  occur  prior  to  attainment  of  the  HEL  if  theoretical  estimate  (4.2)  is 
accurate.  Among  potential  fracture  sites,  failure  on  rhombohedral  planes  is  most 
likely,  but  the  theoretical  limit  from  (4.3)  is  never  attained.  It  thus  appears  likely 
that  nucleated  defects  such  as  intersecting  twins  or  dislocation  pile-ups  could 
produce  stress  concentrations  enabling  fractures  in  the  experiments  (Graham  Sz 
Brooks  1971). 

Collectively,  the  results  for  different  orientations  provide  the  bounds  listed  in 
table  8  (in  the  electronic  supplementary  material).  These  bounds  are  constructed 
from  the  assertion  that  at  least  one  slip  or  twinning  mechanism  occurs  in  a  crystal 
of  each  considered  orientation  upon  attainment  of  the  volumetric  strain 
corresponding  to  the  experimental  HEL  for  that  orientation.  The  lower  bounds 
suggested  in  table  8  for  fracture  on  prism,  pyramidal  and  rhombohedral  planes 
ensure  that  shear  fracture  would  not  occur  prior  to  attainment  of  the  experimental 
HEL  in  any  of  the  considered  crystals.  Bounds  for  fracture  on  the  basal  plane  are 
not  estimated  here.  Owing  to  its  relatively  high  theoretical  and  cohesive  strength 
(table  7),  cleavage  fracture  on  the  basal  plane  seems  unlikely.  Also  shown  are 
bounds  relevant  for  indentation  experiments  at  room  temperature  (Tymiak  Sz 
Gerberich  2007).  The  critical  stress  for  pyramidal  slip  quoted  from  Tymiak  Sz 
Gerberich  (2007)  is  subject  to  large  uncertainty,  since  it  was  obtained  by 
extrapolation,  and  is  treated  here  as  a  lower  bound.  The  bounds  found  via 
consideration  of  shock  physics  experiments  and  nonlinear  elasticity  are  similar  to 
those  from  indentation  (Tymiak  Sz  Gerberich  2007),  apart  from  the  much  higher 
upper  bound  listed  for  basal  slip  by  the  latter  authors.  Confining  pressures  in 
indentation  were  reported  of  the  order  of  100  GPa  (Tymiak  Sz  Gerberich  2007), 
much  higher  than  those  of  the  order  of  10  GPa  considered  here.  Also,  plastic  strain 
rates  lie  in  the  quasi-static  regime  in  indentation,  in  contrast  to  rates  in  excess  of 
105  s-1  attained  in  shock-loading  experiments  (Munson  Sz  Lawrence  1979).  The 
similarity  in  bounds  for  slip  and  twinning  thus  points  to  rate  and  pressure 
independence  of  initial  yield  at  low  to  room  temperatures.  The  pressure 
independence  of  yield  at  confining  pressures  of  the  order  of  1.5  GPa  has  been 
confirmed  experimentally  for  prism  slip  (Castaing  et  al.  1981).  However,  in  shock 
experiments,  plastic  strain  rates  are  not  constant,  but  ramp  up  from  small  values 
at  the  onset  of  yielding  at  the  front  of  the  plastic  wave  to  values  estimated  in  excess 
of  105  s-1  in  the  wake  of  the  plastic  wavefront  (Munson  Sz  Lawrence  1979),  with 
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the  maximum  rate  depending  on  the  material  and  the  severity  of  the  impact 
loading.  While  significant  rate  dependence  of  twinning  has  not  been  observed 
(Scott  &  Orr  1983) ,  minor  rate  sensitivities  of  pyramidal  slip  with  10  ( Tressler  & 

Barber  1974)  and  prism  slip  with  m~  16  (Castaing  et  al.  1981)  have  been  observed 
at  low  loading  rates  and  very  high  temperatures.  Because  strain  rates  are  very 
large  in  impact  experiments,  viscous  drag  acting  on  dislocations  behind  shock 
fronts  (Gilman  1979)  could  probably  elevate  the  flow  stress  for  slip  or  twinning 
such  that  fracture  stresses  on  prism,  pyramidal  and  rhombohedral  planes  would 
be  exceeded  soon  after  initial  yield,  in  the  wake  of  the  plastic  wavefront. 

(c)  Yield 

The  flow  stresses  {5}  entering  (3.45)— (3.47)  are  expressed  as  sums  of  contribu¬ 
tions  of  various  mechanisms,  e.g.  following  Kocks  et  al.  (1975)  or  Clayton  (2005) 

t  =  9s  +  9l ,  9)  =  9jS  +  9jL,  9J  =  9  s  +  9l-  (4-4) 

Here,  gls=g]=  0  reflects  the  initial  yield  stress  for  slip  in  the  parent  crystal  on 
system  i,  and  g\  reflects  long-range  barriers  associated  with  defects  that 
accumulate  during  the  deformation  history.  Analogous  definitions  apply  for 
resistance  to  slip  within  the  twins  and  resistance  to  twinning  in  the  second  and 
third  of  (4.4),  respectively.  Both  terms  in  each  sum  in  (4.4)  depend  on 
temperature;  the  first,  i.e.  initial  yield,  is  addressed  here  for  sapphire,  and  the 
second  term  is  addressed  in  §4 d.  The  initial  yield  stress  depends  on  short-range 
interactions  such  as  strong  Peierls  barriers  in  crystals  with  low  initial  defect 
densities  and  non- metallic  bonds  (Friedel  1964;  Farber  et  al.  1993),  and  at  high 
rates  also  accounts  for  viscous,  phonon  and  electron  drag  (Kocks  et  al.  1975; 
Gilman  1979).  Short-range  barriers  exhibit  a  strong  temperature  dependence 
for  both  slip  and  twinning  in  sapphire  (Lagerlof  et  al.  1994),  reflected 
appropriately  by 

9s  =  9o  exp(— A’0/0m),  gJs  =  ^exp(-A-7 6 / 6 M) ,  (4-5) 

where  gl0  is  an  athermal  yield  stress;  6m  is  the  melting  temperature;  and  A*  is  a 
dimensionless  parameter.  In  the  second  of  (4.5),  analogous  definitions  apply  for 
quantities  associated  with  twin  system  j.  Appropriate  forms  for  the  temperature 
dependence  of  the  second  of  (4.4),  i.e.  for  slip  systems  within  twinned  regions, 
have  not  been  verified  experimentally,  but  it  may  be  reasonable  to  assume  that 
the  short-range  barriers  for  a  slip  system  in  the  parent  are  transferred  to  its 
rotated  counterpart  in  a  twin.  Relation  (4.5)  does  not  conform  to  the  usual 
Arrhenius  form  for  thermally  activated  kinetics  (Kocks  et  al.  1975),  but  can  be 
rationalized  in  terms  of  activation  volumes  for  cross-slip  in  sapphire  (Lagerlof 
et  al.  1994).  Values  of  g0  and  A  are  listed  in  table  9  (in  the  electronic 
supplementary  material)  for  basal,  prism  and  pyramidal  slip  and  basal  and 
rhombohedral  twinning.  Comparisons  with  experimental  data  are  given  in 
figure  5.  Room  temperature  experimental  values  are  taken  from  the  first  column 
of  results  in  table  8,  except  for  pyramidal  slip  whose  strength  is  taken  from  the 
rightmost  column  in  table  8  (Tyrniak  &;  Gerberich  2007).  Pressure  and  rate 
dependence  are  neglected  following  the  discussion  in  §46,  but  could  be  included 
in  a  more  refined  analysis  if  supporting  data  become  available  for  all  slip  and 
twin  system  families.  High-temperature  data  for  basal  and  prism  slip  are 
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0(  K) 

Figure  5.  Experimental  shear  stresses  for  initial  yield  by  slip  or  twinning.  Circles,  basal  slip;  down 
triangles,  prism  slip;  squares,  pyramidal  slip;  diamonds,  basal  twin;  up  triangles,  rlrombohedral 
twin;  solid  line,  model. 

obtained  from  Lagerlof  et  al.  (1994).  Data  for  pyramidal  slip  are  obtained  from 
Tressler  &  Barber  (1974).  For  basal  twinning,  data  are  obtained  from  Castaing 
et  al.  (2004),  and  for  rhombohedral  twinning,  from  Scott  &;  Orr  (1983).  For 
twinning,  the  fits  are  most  appropriate  for  300 <9<  1600  K,  while  for  slip  the 
model  validly  spans  the  range  300<<9<2100  K.  Trend  lines  extended  down  to 
0  K  are  extrapolations.  Predictions  at  300  K  follow 


pyramidal  slip  ^  basal  twin  _  basal  slip  ^  prism  slip  -  rhombohedral  twin 

9s  >55  ~  9s  >  9s  >  9s 


(4.6) 


in  general  agreement  with  the  trends  of  table  8.  A  kink  pair-based  model 
(Rodriguez  et  al.  2008)  may  offer  more  accuracy  for  basal  and  prism  slip  at  high 
temperatures,  albeit  with  more  parameters. 


(d)  Hardening  and  lattice  defects 

Potentially,  the  long-range  barriers  in  each  of  (4.4)  could  depend  in  a  complex 
manner  upon  activity  of  the  slip  or  twin  system  under  consideration  (i.e.  self¬ 
hardening)  as  well  as  the  activity  of  all  other  slip  and  twin  systems  (Christian  & 
Mahajan  1995;  Wu  et  al.  2007).  Considered  here  are  phenomena  for  which 
supporting  data  are  available:  hardening  of  slip  and  hardening  of  twinning  by 
dislocation  accumulation.  To  address  the  former,  let 

n 

9l=  (l-9L/9M)hl^2\yk\,  gl\t=o  =  0,  (4.7) 

k= 1 

where  glM  is  a  saturation  stress  (Wu  et  al.  2007)  and  hl  is  a  hardening  modulus. 
These  vary  with  the  thermodynamic  state  of  the  crystal  as 

9m  =  9m^  p(-«!0/0m),  h*  =  n(hlo-h\d/0M),  (4.8) 
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where  glm,  w1,  hl0  and  h\  are  constants  that  may  differ  among  families  of  slip 
systems.  Relations  (4.7)  and  (4.8)  provide  a  reasonable  fit  to  the  basal  slip  data 
of  Pletka  et  al.  (1977),  with  model  parameters  listed  in  table  9.  Dislocation 
accumulation  takes  place  in  conjunction  with  slip  system  hardening  via 

glL  =  alpb^pT-  pT0,  a1  =  al0-a[d/6M ,  (4.9) 

where  pT 0  is  the  initial  dislocation  density  and  (Xq  and  a\  are  dimensionless 
constants  that  may  vary  among  slip  systems.  Following  Pletka  et  al.  (1977),  (4.9) 
extends  the  theory  of  Taylor  (1934)  of  slip  impedance  from  dislocation 
interactions  on  parallel  planes.  Parameters  are  listed  in  table  9  following  calibra¬ 
tion  to  data  from  experiments  at  a  rate  of  jP  ~js  =  l-33X10_4s_1  (Pletka 
et  al.  1977),  with  yp  the  cumulative  plastic  shear.  Relation  (4.9)  can  be  inverted 
and  used  with  (4.7)  to  provide  an  evolution  equation  for  the  accumulated 
dislocation  line  density  (Clayton  2005).  Note  that  hl  and  glL  are  scaled 
consistently  by  the  rhombohedral  shear  modulus  p  that  depends  on  temperature 
(table  2).  Experiments  also  provide  quantitative  evidence  for  the  hardening  of 
rhombohedral  twin  propagation  by  accumulation  of  forest  dislocations  (Castaing 
et  al.  2002).  Hardening  of  twinning  by  slip  is  captured,  in  general,  by 

n 

9l  =  Y  SlU  =  o  =  0,  (4.10) 

i  =1 

where  h\  is  an  interaction  matrix  relating  the  hardening  rate  of  slip  system  i  to 
that  of  twin  fraction  j.  Parameter  h\  =  hrd\  found  from  a  fit  to  data  of  Castaing 
et  al.  (2002)  is  given  in  table  9,  applicable  only  for  the  hardening  of  rhombohedral 
twinning  by  basal  slip.  Hardening  of  prism  and  pyramidal  slip  systems  and  basal 
twin  systems  by  dislocation  accumulation  has  not  yet  been  addressed.  Twin-twin 
interactions,  hardening  of  slip  systems  within  twins  and  hardening  of  slip  by 
twinning  also  remain  to  be  quantified. 

Stored  elastic  energies  of  dislocation  lines  (energy  per  unit  length)  and  twin 
boundaries  (energy  per  unit  area)  are  written  as 

Wjj  =  ( pb2/(4TrK))ln(R/Rc )  +  W  ~  Kxpb2 /2, 

WT  ~  kFSF/2  =  K2pb/ 2 


where  in  the  first  of  (4.11),  K  accounts  for  the  edge  or  screw  character  of  the 
dislocation  line  as  well  as  elastic  anisotropy  (Eshelby  1949);  R  is  the  radial  distance 
from  the  dislocation  core;  Rc  is  the  cut-off  radius  for  the  dislocation  core;  and  W  is 
a  correction  that  accounts  for  the  core  energy,  line  curvature,  interaction  energies 
from  other  defects  and  boundaries,  and  stacking  faults  associated  with  partial  slip 
dislocations.  As  an  order-of-magnitude  approximation  (Hull  &  Bacon  1984;  Heuer 
et  al.  1998),  here  ki  =  2.0/K  is  assumed.  Elastic  anisotropy  is  neglected  in  this 
context,  and  K  follows  from  (4.1).  Recall  from  §3 d,  and  table  6  that  WgF 
~  1  — 10  J  nr-2  is  a  typical  stacking  fault  energy  associated  with  twin  systems,  from 
which  k2  entering  (3.42)  can  be  quantified  as  Wsp/pb. 

Consider  now  the  net  energy  accumulated  in  the  crystal  during  single  slip  at 
constant  temperature.  From  (4.11),  and  again  assuming  that  rate  effects  on  flow 
stress  are  negligible,  the  cumulative  value  /?  of  the  dissipation  fraction  /?'  in  (3.44) 
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Figure  6.  Model  predictions  for  slip  on  the  basal  plane:  (a)  dislocation  density  accumulation  and 
(6)  cumulative  heat  dissipation  fraction  from  edge  dislocations.  (Solid  line,  0  K;  dotted  line,  300  K; 
small-dashed  line,  600  K;  double  dot-dashed  line,  900  K;  long-dashed  line,  1200  K;  single  dot- 
dashed  line,  1500  K;  medium-dashed  line,  1800  K.) 


can  be  approximated  as 


0  = 


g*dyl  -/xb2 pT/K 


(single  slip). 


(4.12) 


During  basal  slip,  dislocations  are  predicted,  by  inversion  of  (4.9),  to  accumulate  as 
shown  in  figure  6a.  The  corresponding  dissipation  fraction  (3  is  shown  in  figure  6b, 
computed  via  integration  of  (4.7)  and  use  of  (4.9),  (4.12)  and  material  parameters 
obtained  from  fits  to  high-temperature  data  (Pletka  et  al.  1977).  Since  during  the 
experiments  (Pletka  et  al.  1977)  dislocations  of  families  with  various  strengths  and 
orientations  were  generated,  in  the  computations,  the  representative  value  used  for 
b  is  the  rhombohedral  lattice  parameter  and  used  for  g  is  the  rhombohedral  shear 
modulus.  Shown  in  figure  6b  are  results  for  edge  dislocations;  screw  dislocations 
would  result  in  a  decrease  in  energy  storage,  1  —  /?,  by  a  factor  of  K—  0.76.  Note  from 
figure  6  a  that  the  dislocation  densities  increase  with  plastic  shear  from  the  initial 
value  pT0=  10s  m-2  (Pletka  et  al.  1977)  to  saturation  values  in  the  range  1012  m-2 
PT<1014m-2  that  increase  with  decreasing  temperature.  Predicted  values  of  (3 
approach  unity  as  the  temperature  decreases,  primarily  a  result  of  the  high  basal  slip 
resistance  of  sapphire  at  low  temperatures,  which  leads  to  large  dissipation  from 
plastic  work.  At  all  temperatures,  the  stored  energy  appears  small  relative  to  the 
plastic  work,  since  /?>  0.975. 


(e)  Residual  volume  changes 


The  residual  volume  change  J  was  introduced  following  relations  (3.2)  and 
(3.11).  An  estimate  of  this  quantity  from  nonlinear  elasticity  theory  is  (Seeger  &; 
Haasen  1958;  Teodosiu  1982) 


1  fdB 


J  -l  =  Ah/F0=-  --1  4 


B  \  dp 


1  fdG  G 


G  \  dp  B 


(4.13) 


where  G  and  B  are  shear  and  bulk  moduli;  p  is  the  Cauchy  pressure;  and  Ey  and 
Eg  are  dilatational  and  deviatoric  strain  energies  per  unit  volume,  respectively. 
The  volume  change  in  (4.13)  is  measured  between  reference  and  intermediate 
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configurations:  AV=  V  —  V0  with  V  the  volume  in  B.  Partitioning  the  elastic 
energies  per  unit  dislocation  line  length  in  the  first  of  (4.11)  into  dilatational  and 
deviatoric  components  and  neglecting  core  and  interaction  energies,  (4.13)  can  be 
written  as  (Seeger  &  Haasen  1958;  Teodosiu  1982) 


AV 


i 

3 

1 

G 


l-v-2 v2  (dB  \  2-2v  +  2v2 

(1  -v)B  _1  )  +  (1  -v)G 


(screw). 


(edge), 


(4.14) 


Expressions  (4.13)  and  (4.14)  neglect  anisotropy.  An  exact  solution  for  the 
residual  elastic  volume  change,  and  possibly  shape  change,  in  a  highly 
anisotropic  (e.g.  rhombohedral)  lattice  seems  to  require  integration  over  local 
elastic  displacement  gradients  induced  by  defects  within  the  body  (Toupin  Sz 
Rivlin  1960).  Closed- form  expressions  for  residual  elastic  volume  change  exist  for 
anisotropic  single  crystals  of  cubic  symmetry  (Toupin  &  Rivlin  1960),  and  have 
been  applied  in  analysis  of  volume  changes  from  stored  energy  of  cold  work  in 
metal  polycrystals  (Wright  1982).  Here,  following  previous  work  (Seeger  Sz 
Haasen  1958),  isotropic  approximation  (4.14)  is  investigated  as  an  order- 
of-magnitude  estimate,  consistent  with  the  isotropic  approximation  of  dislocation 
energies  in  the  first  of  (4.11).  Table  10  (in  the  electronic  supplementary  material) 
lists  values  of  the  normalized  volume  change  AV/(T62),  where  L  is  the  total 
length  of  dislocation  lines  of  pure  screw  or  pure  edge  character,  and  b  is  the 
Burgers  vector.  For  alumina,  effective  isotropic  constants  (Holm  et  al.  1999)  and 
their  pressure  derivatives  (Sarkar  et  al.  1996)  are  obtained  from  the  literature 
and  are  valid  at  room  temperature.  Values  of  the  normalized  volume  change  in 
AI2O3  from  dislocation  lines  are  positive  in  agreement  with  other  solids  (Seeger  Sz 
Haasen  1958;  Wright  1982),  and  fall  in  between  those  for  metallic  crystal  Cu 
(face-centred  cubic)  and  ionic  crystal  NaCl  (cubic  rock  salt  structure),  with 
properties  of  the  latter  two  solids  given  by  Seeger  Sz  Haasen  (1958)  also  listed  in 
table  10.  Experimental  data  for  inelastic  volume  changes  and  dislocation 
densities  are  available  for  Cu  polycrystals  (Clarebrough  et  al.  1957);  the  range  of 
values  for  the  volume  change  in  the  rightmost  column  of  table  10  corresponds  to 
applied  compressive  strains  ranging  from  0.3-0. 7.  Thus,  the  theory  under¬ 
estimates  the  volume  change  by  a  factor  of  approximately  2-5  for  Cu 
polycrystals.  The  linear  relationship  between  the  dislocation  density  and  volume 
change  in  (4.14)  is  evident.  Wright  (1982)  noted  that  the  ratio  of  dilatational  to 
deviatoric  stored  energy  may  change  during  the  course  of  large  deformation;  this 
could  lead  to  nonlinearity  not  captured  by  either  of  the  estimates  in  (4.14)  alone. 
Volume  changes  predicted  by  (4.14)  remain  small  as  saturation  levels  of 
dislocation  density  implied  by  figure  6a  are  approached.  For  example,  a  density 
of  edge  dislocations  pT—  1014  m-2  would  produce  a  volume  change  A  V/  Vo~4X 
10 _  ,  corresponding  to  the  volume  increase  that  would  result  from  thermal 
expansion  in  alumina  upon  a  temperature  rise  of  approximately  2.5  K.  However, 
according  to  the  theory,  volume  increases  of  the  order  of  1  per  cent  would  be 
achieved  upon  generation  of  immense  dislocation  densities,  e.g.  pT~1016m-2, 
corresponding  to  a  dislocation  spacing  of  the  order  of  10  nm.  Conceivably,  very 
large  dislocation  densities  may  be  required  to  enable  extremely  large  plastic 
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strain  rates  in  shock  physics  experiments  if  the  upper  bound  of  the  dislocation 
velocity  is  limited  to  that  of  transverse  elastic  waves  in  the  crystal  (Kocks  et  al. 
1975).  Generation  of  defect  densities  of  this  magnitude  could  affect  the  measured 
pressure-volume  response;  for  example,  the  pressure  required  to  offset  a  1  per 
cent  volume  increase  would  be  of  the  order  of  0.015=2.5  GPa.  One  should  keep 
in  mind  the  limits  of  the  continuum  elasticity  approximations  used  in  (4.11)  and 
(4.14).  Atomistic  simulations  (Clayton  &  Chung  2006)  may  offer  the  possibility 
of  more  accurate  predictions  of  effects  of  large  defect  densities  on  volume 
changes,  effective  moduli  and  stored  energies. 


5.  Conclusions 

The  theory  developed  here  focuses  on  anisotropic  mechanisms  of  elasticity, 
plastic  slip  and  deformation  twinning.  Large  deformation  theory  is  required  to 
address  finite  shears  arising  from  slip  and  twinning,  lattice  reorientations  arising 
from  twinning  and  nonlinear  elastic  effects.  A  constitutive  framework 
incorporating  internal  state  variables  provides  thermodynamic  relationships 
among  state  variables  and  driving  forces  for  inelastic  deformations.  The 
application  represents  the  first,  to  the  author’s  knowledge,  fully  nonlinear 
anisotropic,  crystal  plasticity-based  model  for  the  ensemble  of  observed 
thermoelastic  and  inelastic  behaviours  of  sapphire.  Inelastic  behaviours  include 
basal,  prism  and  pyramidal  slip  and  basal  and  rhombohedral  twinning.  Estimates 
of  room  temperature  resistances  to  slip  or  twin  initiation  follow  from  nonlinear 
thermoelastic  analysis  and  results  of  prior  shock  physics  experiments  (Barker  &; 
Hollenbach  1970;  Graham  &;  Brooks  1971;  Chen  &  Howitt  1998;  Hare  et  al.  2002; 
Reinhart  et  al.  2006).  These  estimates  are  combined  with  other  existing  data  to 
provide  yield  criteria  over  a  range  of  temperatures.  Hardening  from  dislocation 
accumulation  during  basal  slip  is  quantified,  and  residual  elastic  energies 
associated  with  defects  are  estimated.  The  cumulative  ratio  of  stored  energy  from 
dislocations  to  dissipated  energy  in  basal  slip  is  predicted  at  less  than  3  per  cent 
for  plastic  shears  of  less  than  0.25.  Residual  volume  expansion  from  nucleation  of 
dislocation  lines  is  estimated  using  nonlinear  elasticity  theory  (Seeger  Sz  Haasen 
1958;  Teodosiu  1982).  Such  expansion  is  predicted  to  be  nearly  negligible  in 
sapphire  for  dislocation  line  densities  observed  in  high-temperature  experiments 
(Pletka  et  al.  1977). 
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Table  1  Room  temperature  elastic  moduli 


Parameter 

Value  [GPa] (1) 

Parameter  Value  [GPa]  n> 

Parameter 

Value  [GPa] (1) 

c„ 

498 

Cm 

-3780 

C134 

131 

c12 

163 

Cm 

-1090 

C144 

-302 

c13 

117 

C[  13 

-963 

C155 

-1160 

Cl4 

23 

C114 

-55 

C222 

-4520 

C33 

502 

C123 

-289 

C333 

-3340 

C44 

147 

C[24 

39 

C344 

-1090 

C133 

-922 

C444 

19 

(l)Winey  et  al.  [2001] 

Table  2  Bulk  physical  properties  (atmospheric  pressure) 

Parameter 

Value 

Remarks 

Po 

3980  kg/m3 

room  temperature  mass  density  (1) 

c 

780+0.3(6^-300)  J/kgK 

specific  heat  per  unit  mass  (2) 

«n 

5.0xl0"6  /K 

thermal  expansion  coefficient 1 1  ’ 

«33 

5.7xl0"6  /K 

thermal  expansion  coefficient 1 1  * 

9m 

2325  K 

melting  temperature  (2) 

P 

1 56-70 0I6m  GPa 

rhombohedral  shear  modulus  (3) 

* 1  ’Burghartz  &  Schultz  [1994]  (2)Castanet  [1984]  (3,Zouboulis  &  Grimsditch  [1991] 


Table  3  Slip  systems 


Type 

Direction 

Plane 

Remarks 

Type 

Direction 

Plane  Remarks 

basal 

[1120] 

(0001) 

bidirectional (1) 

pyramidal 

[Tioi] 

(10Tl)  unidirectional <2) 

[12l0] 

(0001) 

[Oil  1] 

(loll) 

[2110] 

(0001) 

[10T1] 

(Oil  1) 

prism 

[llOO] 

(1120) 

bidirectional 

[Tioi] 

(Oil  1) 

[10T0] 

(I2l0) 

[Oil  1] 

(Tioi) 

[01  TO] 

(2110) 

[loll] 

(Tioi) 

(1)Kronberg  [1957];  Snow  &  Heuer  [1973];  Heuer  et  ah  [1998]  <2)Tressler  &  Barber  [1974] 


Table  4  Twin  systems 


Type 

Direction 

Plane 

Shear 

Remarks 

basal 

[T010] 

(0001) 

0.635 

type  II (1) 

[llOO] 

(0001) 

[01  TO] 

(0001) 

rhombohedral 

[Toil] 

(To  12) 

0.202 

type  I (2) 

[01  IT] 

(0ll2) 

[1T0T] 

(1 102) 

(1)Veit  [1921];  Kronberg  [1957]  (2 ’Heuer  [1966] 


Table  5  Slip  dislocations  and  Peierls  stresses 


Plane 

Burgers  vector  (1) 

b  [nm] 

Peierls,  edge  [GPa] 

Peierls,  screw  [GPa] 

basal 

1/3[1 120] 

0.475 

9.0 

16.9 

prism 

1  /  3[1  Too] 

0.274 

0.3 

1.3 

pyramidal 

l/3[Tl01] 

0.512 

0.7 

2.4 

<nKronberg  [1957];  Tressler  &  Barber  [1974];  Lagerlof  et  al.  [1994] 


1 


Table  6  Twinning  partials,  energies,  and  theoretical  nucleation  stresses 


Plane 

Burgers  vector  (1) 

b  [nm] 

Stacking  fault 
energy  [J/m2] <2) 

Nucleation 
stress  [GPa] 

basal 

1/3[I010] 

0.274 

9.0 

7.3 

rhombohedral 

1/ 21.9[1  ToT] 

0.071 

0.7 

3.6 

(1)Kronberg  [1957];  Geipel  et  al.  [1994];  Heuer  et  al.  [1998]  (2,Lagerlof  et  al.  [1984];  Kenway  [1993] 


Table  7  Planar  properties 


Plane 

Full  Burgers 
vector  [nm] ' 1 1 

Spacing  [nm] <2) 

Theoretical  shear 
strength  [GPa] 

Cohesive  energy 

[J/m2] (3) 

basal  {0001} 

0.475 

0.216 

51.5 

>40 

prism  {10  TO} 

0.475 

0.274 

40.5 

7.3 

pyramidal  {10  1 1} 

0.512 

0.392 

30.6 

- 

rhombohedral  {1012} 

0.512 

0.348 

34.4 

6.0 

<nSnow  &  Heuer  [1973]  <2)Geipel  et  al.  [1994];  Lagerlof  et  al  [1994]  <3)Wiederhorn  [1969] 


Table  8  Estimated  bounds  on  shear  strengths  for  slip,  twinning,  and  fracture 


Mechanism 

tc  [GPa]  present  work 

tc  [GPa]  indentation  1 1 1 

basal  slip 

2.7  <  tc<  5.0 

4.8  <  tc<  23.7 

prism  slip 

2.3  <  tc<  5.7 

2.5  <  rc  <  4.1 

pyramidal  slip 

- 

tc  >  7.0 

basal  twin 

tc  >  4.0 

4.8  <  tc<  7.3 

rhombohedral  twin 

1.0  <  tc<  8.1 

1.0  <  rc<7.3 

prism  fracture 

A- >5.9 

~ 

pyramidal  fracture 

Tc  >  6-7 

- 

rhombohedral  fracture 

rc  >  7-9 

- 

<nTymiak  &  Gerberich  [2007] 


Table  9  Slip  and  twinning  parameters 


Parameter 

Value 

Remarks 

Parameter 

Value 

Remarks 

go  [GPa] 

12.7 

basal  slip 

gm  [GPa] 

53.6 

basal  slip 

6.6 

prism  slip 

CO 

10.7 

basal  slip 

10.0 

pyramidal  slip 

ho 

8.0xl0"3 

basal  slip 

11.5 

basal  twin 

h  i 

9.3xl0"3 

basal  slip 

3.4 

rhombohedral  twin 

«o 

0.65 

basal  slip 

A 

8.4 

basal  slip 

«i 

0.70 

basal  slip 

5.6 

prism  slip 

/?,. 

0.42 

rhombohedral  twin 

3.7 

pyramidal  slip 

7.7 

basal  twin 

10.2 

rhombohedral  twin 

Table  10  Effective  isotropic  elastic  properties  and  volume  changes  associated  with  dislocations 


Crystal 

G  [GPa] 

B  [GPa] 

V 

dG/dp 

dB/dp 

AV/(Lb2) 

edge 

AV/(Lb2) 

screw 

AV/(Lb2) 
experiment (3) 

AI2O3 

157 

250 

0.24 

1.7 

4.2 

1.7 

1.0 

— 

NaCl  (2> 

15 

23 

0.24 

2.0 

5.8 

2.5 

1.4 

— 

Cu(2) 

47 

152 

0.36 

0.8 

4.4 

1.0 

0.5 

1. 9-2.6  (edge) 

(1)Sarkar  et  al.  [1996];  Holm  et  al.  [1999]  (2)Seeger  &  Haasen  [1958]  "'Clarebrough  et  al.  [1957] 
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